﻿Imports System.IO
Imports LidarReader
Imports MapWindow.Data

Public Class frmMain

    Dim nNodes As Integer = 5000
    Dim NoiseFactor As Integer = 0.0
    Dim K As Integer = 20 'number of neighbors to search for uphill nodes
    Dim g As New MapWinGIS.Grid
    Dim Outlets As New MapWinGIS.Shapefile
    Dim Nodes() As Node
    Dim cRed As UInt32 = CType(RGB(255, 0, 0), UInt32)
    Dim cBlue As UInt32 = CType(RGB(0, 0, 255), UInt32)
    Dim cGreen As UInt32 = CType(RGB(0, 255, 0), UInt32)
    Dim BasePath As String = "C:\dev\HydroLidar"

    Dim _tree As New MapWindow.Analysis.Topology.KDTree.KDTree(2)



    Private Sub btnLoadDEM_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles btnLoadDEM.Click
        Dim u As New MapWinGIS.Utils
        Dim cs As New MapWinGIS.GridColorScheme
        Dim im As New MapWinGIS.Image
        g.Open(BasePath & "\Data\iron.bgd")
        cs.UsePredefined(g.Minimum, g.Maximum, MapWinGIS.PredefinedColorScheme.SummerMountains)
        im = u.GridToImage(g, cs)
        AxMap1.AddLayer(im, True)

    End Sub

    Private Sub btnLoadOutlets_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles btnLoadOutlets.Click
        Outlets.Open(BasePath & "\Data\outlets.shp")
        Dim hndOutlets As Integer = AxMap1.AddLayer(Outlets, True)
        AxMap1.set_ShapeLayerPointSize(hndOutlets, 10)
        AxMap1.set_ShapeLayerPointColor(hndOutlets, cGreen)
        Dim sf As New MapWinGIS.Shapefile
        sf.Open(BasePath & "\Data\ironw.shp")
        Dim t As Integer = AxMap1.AddLayer(sf, Visible)
        AxMap1.set_ShapeLayerDrawFill(t, False)
    End Sub


#Region "Create Points"
    Private Sub btnCreatePoints_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles btnCreatePoints.Click
        Dim Col As Integer, Row As Integer
        Dim Max As Double = g.Maximum, Min As Double = g.Minimum
        nNodes = nNodes + Outlets.NumShapes
        AxMap1.NewDrawing(MapWinGIS.tkDrawReferenceList.dlSpatiallyReferencedList)

        ReDim Nodes(nNodes - 1)

        For i = 0 To nNodes - 1
            Nodes(i) = New Node
            If i < Outlets.NumShapes Then       'use specific outlet point
                Nodes(i).X = Outlets.Shape(i).Point(0).x
                Nodes(i).Y = Outlets.Shape(i).Point(0).y
            Else                                'use random point
                Nodes(i).X = RandBetween(g.Header.XllCenter, g.Header.XllCenter + g.Header.NumberCols * g.Header.dX)
                Nodes(i).Y = RandBetween(g.Header.YllCenter, g.Header.YllCenter + g.Header.NumberRows * g.Header.dY)
            End If
            Dim dbl(1) As Double
            dbl(0) = Nodes(i).X
            dbl(1) = Nodes(i).Y
            _tree.Insert(dbl, Nodes(i))
            g.ProjToCell(Nodes(i).X, Nodes(i).Y, Col, Row)
            Nodes(i).Z = AddNoise(g.Value(Col, Row), Max - Min)
            AxMap1.DrawPoint(Nodes(i).X, Nodes(i).Y, 3, cRed)
            'AxMap1.AddLabel(0, Nodes(i).Z, cRed, Nodes(i).X, Nodes(i).Y, MapWinGIS.tkHJustification.hjCenter)
        Next

    End Sub

    Private Function RandBetween(ByVal Low As Double, ByVal High As Double) As Double
        Return (High - Low) * Rnd() + Low
    End Function

    Private Function AddNoise(ByVal Value As Double, ByVal Range As Double) As Double
        Dim Noise As Double
        Noise = Rnd() * Range / 100 * NoiseFactor
        Noise -= Noise / 2
        Return Value + Noise
    End Function

#End Region



    Private Sub btnTry_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles btnTry.Click
        'use "-" to indicate outlets at virtual points
        Dim i As Integer
        For i = 0 To 4
            ProcessNode(Nodes(i))
        Next

        'instead use some seed nodes, all the node in the Nodes will be check if they have DownNodes and the DownNodes should 
        'obey some rules

        'For i = 0 To Nodes.Length - 1
        '    GetDownNode(Nodes(i))
        'Next
    End Sub
    Private Sub GetDownNode(ByVal CurrentNode As Node)
        Dim Element(1) As Double
        Dim slope As Double, i As Integer
        Element(0) = CurrentNode.X
        Element(1) = CurrentNode.Y
        Dim Neighbors() As Object
        Neighbors = _tree.Nearest(Element, K)
        For i = 0 To K - 1
            Dim Neighbor As Node
            Neighbor = CType(Neighbors(i), Node)
            If Neighbor.Z > CurrentNode.Z Then
                'compute slope to this uphill node and see if it is steeper than any existing slope stored at that node
                Dim oldSlope As Double = Neighbor.Slope()
                slope = Neighbor.GetSlope(CurrentNode)
                If slope > oldSlope Then
                    Neighbor.DownNode = CurrentNode
                End If
            End If
        Next
        'CurrentNode.Visited = True
        'draw this node on the map
        AxMap1.DrawCircle(CurrentNode.X, CurrentNode.Y, 4, cRed, False)
    End Sub
    Private Sub ProcessNode(ByVal CurrentNode As Node)
        'compute distance to every node
        Dim d(nNodes - 1) As Double, keys(nNodes - 1) As Integer
        Dim slope As Double, i As Integer
        Dim Location(1) As Double
        Location(0) = CurrentNode.X
        Location(1) = CurrentNode.Y
        Dim Neighbors() As Object
        Neighbors = _tree.Nearest(Location, K)

        'check nearest K nodes to see if any are higher in elevation than the current node.
        For i = 0 To K - 1
            Dim Neighbor As Node
            Neighbor = CType(Neighbors(i), Node)

            If Neighbor.Z > CurrentNode.Z Then
                'compute slope to this uphill node and see if it is steeper than any existing slope stored at that node
                Dim oldSlope As Double = Neighbor.Slope()
                slope = Neighbor.GetSlope(CurrentNode)
                If slope > oldSlope Then
                    Neighbor.DownNode = CurrentNode
                End If
                'now recursively do this uphill node... 
                If Neighbor.Visited = False Then
                    ProcessNode(Neighbor)
                End If
            End If
        Next
        CurrentNode.Visited = True
        'draw this node on the map
        AxMap1.DrawCircle(CurrentNode.X, CurrentNode.Y, 4, cRed, False)
    End Sub

    Private Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click
        'draw lines
        Dim i As Integer
        For i = 0 To nNodes - 1
            If Not Nodes(i).DownNode Is Nothing Then
                AxMap1.DrawLine(Nodes(i).X, Nodes(i).Y, Nodes(i).DownNode.X, Nodes(i).DownNode.Y, 1, cBlue)
            End If
        Next
        '        AxMap1.CursorMode = MapWinGIS.tkCursorMode.cmNone
        '       AxMap1.SendMouseDown = True
    End Sub

    Private Sub AxMap1_MouseDownEvent(ByVal sender As System.Object, ByVal e As AxMapWinGIS._DMapEvents_MouseDownEvent) Handles AxMap1.MouseDownEvent
        Dim xy(1) As Double
        Dim px As Double, py As Double
        Dim n() As Object, nn As Node
        Dim i As Integer
        AxMap1.PixelToProj(e.x, e.y, px, py)
        xy(0) = px
        xy(1) = py
        n = _tree.Nearest(xy, 10)
        For i = 0 To n.GetUpperBound(0)
            nn = CType(n(i), Node)
            AxMap1.DrawCircle(nn.X, nn.Y, 5, cBlue, True)
        Next i

    End Sub

    Private Sub frmMain_Load(ByVal sender As Object, ByVal e As System.EventArgs) Handles Me.Load
        AxMap1.NewDrawing(MapWinGIS.tkDrawReferenceList.dlSpatiallyReferencedList)
    End Sub

    Private Sub ReadLAS_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ReadLAS.Click
        Dim filename As String = BasePath & "\Data\1.las"
        Dim LidarStream As New FileStream(filename, FileMode.Open, FileAccess.Read, FileShare.Read)
        Dim br As New BinaryReader(LidarStream)
        Dim hdr As New LasHeader
        hdr.Read(br)
        br.BaseStream.Seek(hdr.PointDataOffset, SeekOrigin.Begin)
        Dim LidarFeatureSet As New MapWindow.Data.FeatureSet

        Dim i As Integer
        nNodes = hdr.Count
        ReDim Nodes(nNodes - 1)
        For i = 0 To hdr.Count - 1
            Dim lidpnt As New LidarPoint
            lidpnt.ReadPoint(br, hdr)
            Dim pntLidar As New MapWindow.Geometries.Point
            pntLidar.X = lidpnt.X
            pntLidar.Y = lidpnt.Y
            pntLidar.Z = lidpnt.Z
            Dim f As New MapWindow.Data.Feature(pntLidar)
            LidarFeatureSet.AddFeature(f)
            Nodes(i) = New Node
            Nodes(i).X = pntLidar.X
            Nodes(i).Y = pntLidar.Y
            Nodes(i).Z = pntLidar.Z
            Dim dbl(1) As Double
            dbl(0) = Nodes(i).X
            dbl(1) = Nodes(i).Y
            _tree.Insert(dbl, Nodes(i))
        Next
        LidarFeatureSet.SaveAs(BasePath & "\Data\TestLidar.shp", True)
        Dim SL As New MapWinGIS.Shapefile
        SL.Open(BasePath & "\Data\TestLidar.shp")
        AxMap1.AddLayer(SL, Visible)


       

    End Sub
End Class
